;+
; NAME:
;       MARGINALIZE_2D_PARAMS
;
; PURPOSE:
;       Find the mariginalized distributions for 2-D joint-confidence
;       parameter sets (from chi-squared analyses).
;
; CATEGORY:
;       Statistics, or something
;
; CALLING SEQUENCE:
;       result = MARGINAZLIZE_2D_PARAMS, A, B, CHISQ
;
; INPUTS:
;       A      -- The array of A parameters
;       B      -- The array of B parameters
;       CHISQ  -- The array of chisq values
;
; OPTIONAL INPUTS:
;       NONE
;
; KEYWORD PARAMETERS:
;       NONE
;
; OUTPUTS:
;       RESULT -- Structure containing the following elements:
;                 .A_value -- The best value for A (marginalized over B)
;                 .A_minus -- The negaitve "1-sigma" error bar for A
;                 .A_plus  -- The positive "1-sigma" error bar for A
;                 .B_value -- The best value for B (marginalized over A)
;                 .B_minus -- The negaitve "1-sigma" error bar for B
;                 .B_plus  -- The positive "1-sigma" error bar for B
;
; OPTIONAL OUTPUTS:
;       NONE
;
; MODIFICATION HISTORY:
;
;       Created:  04/21/11, TPEB -- Initial version.
;       Modified: 05/02/11, TPEB -- Changed the internal workings of
;                                   the routine to deal with
;                                   delta_chisq so that likelihoods
;                                   are not machine 0.
;       Modified: 05/18/11, TPEB -- Changed uncertainties in A & B
;                                   parameters to be interpolated
;                                   values rather than grid values.
;
;-

FUNCTION MARGINALIZE_2D_PARAMS, A, B, chi, logy
  
  ;; Create new variable so we don't &*^$@* up the original
  chisq = chi
  
  ;; Create the output structure
  result = CREATE_STRUCT('A_value',0.,'A_minus',0.,'A_plus',0.,$
                         'B_value',0.,'B_minus',0.,'B_plus',0.)
  
  ;; First make chisq FINITE!
  find = WHERE(FINITE(chisq))
  A     = A[find]
  B     = B[find]
  chisq = double(chisq[find])
  
  ;; Next, work exclusively in DELTA_CHISQ (else probabilities may get
  ;; too small to work with)
  chisq = chisq - MIN(chisq)
    
  ;; Convert chisq into likelihoods!
  ll = exp(-chisq/2.d)
    
  ;; Identify the VALUES of the A & B parameters, and get sizes
  A_vals = A[UNIQ(A, SORT(A))]
  B_vals = B[UNIQ(B, SORT(B))]
  n_A = n_elements(A_vals)
  n_B = n_elements(B_vals)
  
  ;; First marginalize over A...
  A_ll = dblarr(n_A)
  
  FOR i=0L, n_A-1 DO BEGIN
     ind = WHERE(A EQ A_vals[i], nind)
     IF nind GT 0 THEN $
        A_ll[i] = TOTAL(ll[ind],/DOUBLE,/NAN)
  ENDFOR

  ;; Reconvert to a chisq (--> redch = delta chi^2)
  A_chisq = -2.d * ALOG(A_ll / TOTAL(A_ll,/DOUBLE,/NAN))
  redch = A_chisq - MIN(A_chisq, minind)
  
  result.A_value = median(A_vals[WHERE(redch EQ 0.)])

  ;; Interpolate to find where redch = 1
  lin = WHERE(A_vals LE A_vals[minind],nlin)
  IF nlin LE 1 THEN result.A_minus = 0. ELSE $
     result.A_minus = ABS( (INTERPOL(A_vals[lin],redch[lin],1.) > MIN(A_vals))$
                           - result.A_value)
  
  hin = WHERE(A_vals GE A_vals[minind],nhin)
  IF nhin LE 1 THEN result_A_plus = 0.  ELSE $
     result.A_plus  = ABS( (INTERPOL(A_vals[hin],redch[hin],1.) < MAX(A_vals))$
                           - result.A_value)
  
  plot,A_vals,A_ll,color=120,thick=5,linestyle=0,xst=5,yst=4  
  vline,A_vals[WHERE(A_ll EQ MAX(A_ll))],thick=5,linestyle=4,color=120  
  
  
  ;; next marginalize over B...
  B_ll = dblarr(n_B)
  
  FOR i=0L, n_B-1 DO BEGIN
     ind = WHERE(B EQ B_vals[i], nind)
     IF nind GT 0 THEN $
        B_ll[i] = TOTAL(ll[ind],/DOUBLE,/NAN)
  ENDFOR
  
  ;; Reconvert to a chisq (--> redch = delta chi^2)
  B_chisq = -2.d * ALOG(B_ll / TOTAL(B_ll,/DOUBLE,/NAN))
  redch = B_chisq - MIN(B_chisq, minind)
  
  result.B_value = B_vals[WHERE(redch EQ 0.)]

  ;; Interpolate to find where redch = 1
  lin = WHERE(B_vals LE B_vals[minind],nlin)
  IF nlin LE 1 THEN result.B_minus = 0. ELSE $
     result.B_minus = ABS( (INTERPOL(B_vals[lin],redch[lin],1.) > MIN(B_vals))$
                           - result.B_value)
  
  hin = WHERE(B_vals GE B_vals[minind],nhin)
  IF nhin LE 1 THEN result_B_plus = 0.  ELSE $
     result.B_plus  = ABS( (INTERPOL(B_vals[hin],redch[hin],1.) < MAX(B_vals))$
                           - result.B_value)
  
  
  plot,B_ll,B_vals,color=220,thick=5,linestyle=0,xst=4,yst=5,ylog=logy
  vline,/horiz,B_vals[WHERE(B_ll EQ MAX(B_ll))],thick=5,linestyle=4,$
                      color=220
  
  
  
  RETURN,result
END
